function
   Nuc0(ξi, n) v = CIJ[n, ξi,:]
   r = zeros(1,10*ndp)
   r[1 : 10 : end] = v
   return r
end

function
   Nθc(ξi, n) v = CIJ[n, ξi,:];
   r = zeros(1,10*ndp)
   r[2 : 10 : end] = v
   return r
end

function
   Nwc(ξi, n) v = CIJ[n, ξi,:]
   r = zeros(1,10*ndp)
   r[3 : 10 : end] = v
   return r
end

function Nus0(ξi, n)
   v = CIJ[n, ξi,:]
   r = zeros(1,10*ndp)
   r[4 : 10 : end] = v
   return r
end

function Nus1(ξi, n)
   v = CIJ[n, ξi,:]
   r = zeros(1,10*ndp)
   r[5 : 10 : end] = v
   return r
end

function Nus2(ξi, n)
   v = CIJ[n, ξi,:]
   r = zeros(1,10*ndp)
   r[6 : 10 : end] = v
   return r
end

function Nus3(ξi, n)
   v = CIJ[n, ξi,:]
   r = zeros(1,10*ndp)
   r[7 : 10 : end] = v
   return r
end

function Nws0(ξi, n)
   v = CIJ[n, ξi,:]
   r = zeros(1,10*ndp)
   r[8 : 10 : end] = v
   return r
end

function Nws1(ξi, n)
   v = CIJ[n, ξi,:]
   r = zeros(1,10*ndp)
   r[9 : 10 : end] = v
   return r
end

function Nws2(ξi, n)
   v = CIJ[n, ξi,:]
   r = zeros(1,10*ndp)
   r[10 : 10 : end] = v
   return r
end

Nuc(ξi, z) = Nuc0(ξi, 0) - z*Nθc(ξi, 0)
Nus(ξi, z) = Nus0(ξi, 0) + z*Nus1(ξi, 0) + z^2*Nus2(ξi, 0) + z^3*Nus3(ξi, 0)
Nws(ξi, z) = Nws0(ξi, 0) + z*Nws1(ξi, 0) + z^2*Nws2(ξi, 0)

Nucd(ξi, z) = Nuc0(ξi, 1) - z*Nθc(ξi, 1)
Nusd(ξi, z) = Nus0(ξi, 1) + z*Nus1(ξi, 1) + z^2*Nus2(ξi, 1) + z^3*Nus3(ξi, 1)
Nwsd(ξi, z) = Nws0(ξi, 1) + z*Nws1(ξi, 1) + z^2*Nws2(ξi, 1)

function Nϵsx(xi, xj, ξi, z)
   invj = 2/(xj - xi)
   return invj*Nusd(ξi, z)
end

function Nϵsz(ξi, z)
   return Nws1(ξi, 0) + 2*z*Nws2(ξi, 0)
end

function Nγs(xi, xj, ξi, z)
   invj = 2/(xj - xi)
   return Nus1(ξi, 0) + 2*z*Nus2(ξi, 0) + 3*z^2*Nus3(ξi, 0) + invj*Nwsd(ξi, z)
end

function Nϵc(xi, xj, ξi, z)
   invj = 2/(xj - xi)
   return invj*Nucd(ξi, z)
end

function Nγc(xi, xj, ξi)
   invj = 2/(xj - xi)
   return invj*Nwc(ξi, 1) - Nθc(ξi, 0)
end

function Nucs(ξi, h1, h2)
   return Nuc(ξi, -h1) - Nus(ξi, h2)
end

function Nwcs(ξi, h2)
   Nwc(ξi, 0) - Nws(ξi, h2)
end

function BMatrixC(xi, xj, ξi, z)
   return [Nϵc(xi, xj, ξi, z); Nγc(xi, xj, ξi)]
end

function BMatrixS(xi, xj, ξi, z)
   return [Nϵsx(xi, xj, ξi, z); Nϵsz(ξi, z); Nγs(xi, xj, ξi, z)]
end
